unit Unit1;

interface

type
 TPDFFunction = function(a,x: double): double;
 TCDFFunction = TPDFFunction;
 function NeumanRand: double; export;
 function InverseRand: double; export;
 procedure InitRandomFunctions(pfun: TPDFFunction; cFun: TCDFFunction);export;
 procedure SetParameters(paramA, right, maxPDF: double; valuesCount: integer);export;

 var
  pdfFun: TPDFFunction;
  cdfFun: TCDFFunction;
  paramA_: double;
  valuesCount_: integer;
  right_:double;
  maxPDF_:double;
  pdfVaues: array of double;
  xValues: array of double;
implementation

 procedure SetParameters(paramA, right, maxPDF: double; valuesCount: integer);
 var intervalLen: double;
     i,j, steps: integer;
 begin

  paramA_:=paramA;
  valuesCount_:=valuesCount;
  right_:=right;
  maxPDF_:=right;
  SetLength(pdfVaues, valuesCount_);
  SetLength(xValues, valuesCount_);

  intervalLen:=right_/valuesCount_;
  for i:= 0 to valuesCount_-1 do
   begin
    pdfVaues[i]:=cdfFun(paramA_, j*intervalLen);
{    for j:=0 to i-1 do
     pdfVaues[i]:=pdfVaues[i] + pdfFun(paramA_, j*intervalLen)*intervalLen;
     xValues[i]:=i*intervalLen;
   end;
}
 end;

 function NeumanRand: double;
 var x,y:double;
 begin
  result:=0;
   while (true) do
   begin
    y:=maxPDF_*random;
    x:=right_*random;
    if y < pdfFun(paramA_, x) then
     begin
      result:=x;
      break;
     end;
   end;
 end;

 function InverseRand: double;
 var i: integer;
     y: double;
 begin
  y:=random;
  for i:= 0 to valuesCount_-1 do
  begin
   if (y - pdfVaues[i]) < 0 then
    begin
     result:=xValues[i];
     break;
    end; 
  end;
 end;

procedure InitRandomFunctions(pfun: TPDFFunction; cFun: TCDFFunction);
begin
 pdfFun:= pfun;
 cdfFun:= cfun;
end;

end.
